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Abstract. For a next-nearest neighbour pair interaction model in a periodic 
domain, a priori and a posteriori analyses of the quasinonlocal quasicontinuum 
method (QNL-QC) are presented. The results are valid for large deformations 
and essentially guarantee a one-to-one correspondence between atomistic solutions 
and QNL-QC solutions. The analysis is based on truncation error and residual 
estimates in negative norms and novel a priori and a posteriori stability estimates. 
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O 1. Introduction 
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Quasicontinuum (QC) methods [TOj [T7J HU [20j 126] are a class of prototypical 
schemes for coupling atomistic models of solids to continuum elasticity models. They 
combine, in principle, the accuracy of atomistic models for defects with the efficiency 
of continuum models for long-range elastic effects. For further detail on the historical 
development of the QC method see [171 US] . 

Not all QC methods are equally reliable in modelling large atomistic systems. 
Recent analyses of QC methods [SI El El El [91 H21 [J6l H9] have identified various 
sources of inconsistencies (in the sense of numerical analysis; see also Section [3]) and 
instabilities due to inadequate treatments of the interface. The two methods that 
have emerged as the most promising candidates are the force-based QC method [I] 
^ and the quasinonlocal QC (QNL-QC) method [26] (and its generalization [TU] to 

C""-. longer interaction ranges). 

All contributions to the analysis of QC methods cited above consider only lin- 
earized problems, except for [19] where small deformations from a uniform reference 
state are admitted. The main purpose of the present paper is to introduce a simple 
Q\ yet powerful analytical framework that allows a generalization of these results to 

genuinely nonlinear situations. To present the new ideas that are required to achieve 
this in the simplest possible setting, we will focus on a one-dimensional model prob- 
^ lem where the atomistic model is formulated in terms of a second neighbour pair 

interaction energy with applied dead loads. For this model, the QNL-QC method of 
Shimokawa et al. [26J was previously analyzed in [HI [HJ] . However, the techniques 
used there would be difficult to extend to the large deformation regime, for which an 
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entirely new approach will be presented here, based on consistency error estimates in 
negative Sobolev norms and sharp a priori and a posteriori stability estimates. More- 
over, care will be taken throughout this paper that the techniques can, in principle, 
be extended to higher dimensions and to situations with defects. This will primar- 
ily be achieved by modifying ideas from [22J (where a Galerkin projection without 
continuum approximation is analyzed), so that no higher regularity results for the 
hessian operator of the atomistic or QC energy functionals are required. 

The present paper also develops the first residual-based a posteriori error analysis 
for the QNL-QC method. Theorem [9] is the first a posteriori error estimate for 
the QNL-QC method which rigorously establishes the existence of exact solutions 
for which the estimate applies. The application of goal oriented a posteriori error 
control to QC methods was pioneered in (TJ [2j 125] . 

At several places in the paper interesting and challenging open problems are dis- 
cussed. 

Finally, it should be remarked that the purpose of this paper is to present a new 
framework for the analysis of the QC method in the simplest possible setting that 
still allows the main features one observes in applications: non-smooth solutions and 
large deformations. Numerical experiments illustrating the results presented here, 
as well as further interesting observations will be presented after introducing finite 
element coarse graining into the framework [23] . 



Outline. We begin, in Section |2.1[ by presenting an atomistic model problem in a 



periodic formulation that avoids the difficulties posed by boundaries. In Section 2.2 
we formulate the QNL-QC approximation without coarse graining. For an analysis 
of the QNL-QC including finite element coarse graining of the continuum region see 
|23j . In the remainder of section [2] we discuss some features of the atomistic model 
and introduce the necessary technical background for the subsequent analysis. 

In Section [3] we discuss the concepts of truncation error and residual in the context 
of the QC method. Previous work [5] has analyzed the truncation error in weighted t v - 
norms, which leads to suboptimal truncation error estimates. As a result, somewhat 
subtle and technically demanding ideas were required to still obtain (quasi-)optimal 
error estimates. By contrast, we will see here that if the "correct" norm is chosen, 
namely a negative Sobolev norm, then these complications can be completely avoided. 
In [19], an interesting variant of a negative norm was used that is useful for the 
nonlinear analysis but does not lead to (quasi-) optimal error estimates (see Remarks 
gjandg). 

The subject of Section [4] is a careful and general stability analysis of the QNL-QC 
method. First, an idea from [9] will be extended to obtain a priori stability results 
for large deformations. The main novel contribution in this section, however, is the a 
posteriori stability result in Theorem [6] This section concludes with a brief discussion 
of the difficulties encountered when defects are present, and how the results might 
be generalized. 

Finally, in Section [5] we combine the consistency and stability analyses of the 
previous sections to obtain a priori as well as a posteriori existence results and error 
estimates. The techniques used here are not too dissimilar from [19], however, the 
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point of view taken is a very different one and leads to a different set of existence 
results. Moreover, the sharper consistency and stability results of the present paper 
lead to error estimates that are truly (quasi-)optimal in the atomistic spacing and 
the smoothness of the atomistic solution. See Remarks |9] and [10] for more detail. 

2. Atomistic Model and QNL-QC Approximation 

2.1. An atomistic model problem. In atomistic models surfaces create boundary 
layers in the deformation field and must therefore be considered defects in the lattice. 
We avoid this issue, by working in an infinite lattice eZ, where e > is the reference 
lattice spacing. Since the functional analysis becomes unnecessarily difficult in this 
infinite domain we will, moreover, restrict the set of admissible deformations to those 
which are iV-periodic displacements from the reference lattice, that is, we define 

& = Fx + W, 

where = e(, for £ G Z, and F > is a macroscopic deformation gradient, and 
where 

^ = {m6K z : u i+N = m 5 V£ G Z, and J2f=i u H = °}- 

The condition Ylf=i = is an artefact of the periodic boundary conditions, and 
ensures locally unique solvabililty of the equilibrium equations to the models that 
are introduced below. We set e = l/N throughout so that the reference length of 
the period is one. 

We denote the algebraic dual of by . Equipped with the weighted £ 2 -inner 
product 

N 

(v,w) =e^2v/:w^ 
€=i 

% becomes a Hilbert space, and hence, any T G can be represented by a function 
z T G % via T[v] = (z T ,v) for v G % . Thus, even though forces are usually best 
understood as elements of we will usually identify them with "functions" / G ^/ . 

The stored energy per period of a deformation y G W is given by a next-nearest 
neighbour pair interaction model, 

N 

*(*):= + ^ + (1) 

where v'^ = — t^-i) for v G M z , and where is a Lennard- Jones type interaction 

potential which satisfies the following properties: 

(01) <p G C 3 ((0,+oo];M), and 

(02) there exists r* > such that is convex in (0, r*) and concave in (r*, +oo). 

By G C 3 ((0, +oo]; K) we mean that and its first three derivatives are bounded 
in any interval {5, +oo), 5 > 0. Implicitly, we actually assume that 0(r) and its 
derivatives decay rapidly to zero as r — > oo. This justifies the next-nearest neighbour 
interaction model. 
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We assume, for simplicity, that the atomistic system is subjected to a dead load 
/ G so that the total energy of a deformation y G & is given by 

&°\y) = *(i/) -</,y>. 

Our goal is to solve the minimization problem 

y G argmin$ tot (^), (2) 

where argmin denotes the set of local minimizers. The first order necessary optimality 
condition (or, Euler-Lagrange equation) for (|2]) is 

D^ tot (y)[u] = VueW. (3) 

A deformation y G W satisfying ^ is called an equilibrium (or critical point). If an 
equilibrium y also satisfies the sufficient second order optimality condition 

D 2 <5> tot {y)[u,u] > V«G^\ {0}, (4) 

then we say that y is a strongly stable equilibrium. A strongly stable equilibrium is 
an isolated local minimizer of $ tot in W . 

2.2. The quasi-nonlocal QC method. If a deformation y G IV is smooth, in the 
sense that ~ for all £, then we can approximate the second-neighbour terms 
by 

*(l4+l4fi) ~|{0(2^) + 0(2^+i)}> (5) 
which leads to an approximate stored energy functional 

N N 
$ cb (y) = 8 {<f>(vd + 0(V € )} = ^0cb(^), 

where c b(^) = + <p(2r) is called the Cauchy-Born stored energy function. 
The Cauchy-Born stored energy functional $ c b is local and its minimization can 
be achieved efficiently by means of finite element methodology. (Note, though, that 
in ID this is not really an issue |22j.) Moreover, it is considered an excellent ap- 
proximation to the atomistic model in a smooth deformation regime (see [TTJ for a 
rigorous theory of the small deformation regime). 

If a minimizer of the atomistic model, which we are trying to compute, has defects 
in a small localized region but is smooth elsewhere then we should couple the atom- 
istic model to the continuum Cauchy-Born model. The QNL-QC method [22] is the 
simplest energy-based coupling method that is consistent for next-nearest neighbour 
models (see Section [3] for a discussion of consistency). For the simple pair interaction 
energy we use in this paper it can be defined as follows. 

Let A C {1, . . . , N} be the atomistic region and let C = {1, . . . , N} \ A be the 
continuum region. The nearest neigbour terms are left unchanged, while, for £ G A, 
the next-nearest neighbour interaction <j)(y^ + y'g+i) is approximated as in Q. This 
procedure leads to the QNL-QC stored energy functional 

TV 

$ qc (y) = e + £ E Me + »W + £ E + We+i) } • (6) 

e=l £eA tec 
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Figure 1. Visualization of the atomistic region, the continuum re- 
gion, the interface sets, as well as their counterparts for bounds, defined 
in Section 12^21 

Upon defining the total energy 

in the QNL-QC method, we solve the minimization problem 

Vqc G argmin $*° t (^). (7) 
The associated first order optimality condition in variational form is 

£>$£(y qc )[u] = Vwer (8) 

Remark 1. For practical reasons, it is important that this approximation can be 
rewritten in the form 

**(y) = ^(y) + *£(y) + ^(y), 

where <3>q C , and $q C are the contributions from the atomistic region, the interfa- 
cial region, and the continuum region, respectively, and so that is precisely the 
atomistic model restricted to (a subset of) A, that $q C is precisely the Cauchy-Born 
approximation restricted to (a subset of) C, and $q C is an interface correction that 
is cheap to compute. This formulation allows for finite element coarsening of the 
continuum region. For more details on this splitting and its importance, see [23j. For 
the purposes of the analysis presented in this paper, however, ^ is the most useful 
formulation. □ 
For future reference, we define the left- and right-interface sets (left and right are 
understood in relation to the atomistic region(s)) as follows, 

Xi= {£eC:£ + l eA#] and l r = {£ G C : f - 1 G A#], (9) 

where A# denotes the periodically extended atomistic region, 

A # = {J{Nj+A). 

The union I = If U I r is simply called the interface. A visualization of these, and 
the following, definitions is shown in Figure [Tj 
We assume throughout that 

T £ nJ r = 0. (10) 

This means, that any connected component of the continuum region must contain at 
least two atoms. 
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In addition, since we will often write expressions in terms of bonds rather than 
atoms, it is sometimes more convenient to use the following variants of these sets: 

j;=J £ + l, j;=X r , A' = A\l' e , C' = CUl' e , and l' = l' r Ul' e . 

These definitions are made mostly for the sake of an intuitive notation in the main 
results. On a first reading, it may be best to ignore these subtle differences and 
simply bear in mind that the sets A', etc., are variants of the originals that allow for 
an attractive presentation. 

2.3. A bound on next-nearest neighbour interactions. It will be a crucial in- 
gredient in the analysis of the QNL-QC method to assume that next-nearest neigh- 
bour interactions remain in the concave region of the interaction potential, that is, 
we assume that the solutions of the atomistic model and of the QNL-QC method 
satisfy 

y[>rj2 V£e{l,...,iV}. (ll) 

For most interesting interaction potentials, the case < r*/2 can only be achieved 
under extreme compressive forces that will not usually be observed in experiments. 
As a matter of fact, one can safely assume that under such extreme conditions a 
classical potential (as opposed to a potential based on quantum mechanics) and zero 
temperature statics constitute an inappropriate model to begin with. 

Note, moreover, that the condition y^ > r*/2 also ensures that and its derivatives 
satisfy uniform bounds. For future reference, we define the constants 

Cj(s) := sup0^(r) for s > and for k — 0, 1, ... , 

r>s 

where <j>^ denotes the j'th derivative of 4>. 

2.4. Further notation and an auxiliary result. The second and third finite dif- 
ferences are defined, for t; 6 l z , by 

v'l := e~ 2 (v^ + i - 2v£ + and 
v'l' := e^ 3 (v H i - 3v ( + 3w ? _i - v^ 2 )- 

For «6f , the functions v', v", and v'" are understood as iV-periodic functions. We 
normally associate the function value v% and the second difference v'^ with the atom 
£, while the first difference v'^ and the third difference v'^' are normally associated 
with the cell, or bond, (£ — 1,£). 

In addition to the weighted ^-inner product we define the weighted £ p -norms 

„ „ / (sEJLiW) 1 ", i<p<oo, 

\\ v \\e p '■= \ ' ii 

" " s I max \ve\, p = oo, 

We normally drop the subscript in if and write £°° instead.. 

The space will often be equipped with the Sobolev-like norms 

||w||^/i, P = H^'ll^f, for t)£f and p G [1, oo]. 
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The space % equipped with the < ^' 1 ' p -norm is then simply denoted % l,v . For p' = 
p/(p — 1), the norm on the topological dual '%f~ l,p := (^ 1,p )* of '%f l,p is denoted 

||T||^-i, P := ||T|| ( ^i, y) « = sup T[v], for T G ^~ llP . 

IMI<W =1 

If \1/ : W — > R is (Frechet) differentiable then its first variation is denoted -D\l/ and 
is understood as a nonlinear map from W to that is, for y G D^>(y) G 
f i — > D^f(y)[v]. Similarly, the second variation D 2 ^! is understood as a nonlinear map 
from W to L(W,W*), that is, for |/ef, £> 2 ^(?/) G L(^, ^*). Equivalently, L> 2 #(?/) 
can be understood as a symmetric bilinear form on % x (i>, w) i— > D 2 ^(y)[i;, itj]. 
We will henceforth make no distinction between these different interpretations and 
use whichever is most convenient at any given moment. 

Lemma 1 (Inverse Function Theorem). Let X, Y be Banach spaces, A an open 
subset of X , and let & : A —>■ Y be Frechet differentiable. Suppose that xq G A 
satisfies the conditions 

||^Oo)||y < v, \\D&(x a y l \\ L( Y,x) < o-, 

B x (x ,2r)a) C A, 

\\D&(xi) - D^(x 2 )\\l(x,y) < L\\xi - x 2 \\x for \\xj - x \\ x < 2rja, 
and 2La 2 i] < 1, 

then there exists x G X such that ^(x) = and \\x — Xq\\x — 2r\a . 

Proof. The result follows, for example, by applying Theorem 2.1 in [21] with the 
choices R = 2r/a, u(xq,R) = LR and lj(xq,R) = \LR 2 . Similar results can be 
obtained by tracking the constants in most proofs of the inverse function theorem, 
and assuming local Lipschitz continuity of . □ 

3. Consistency 

It was pointed out in [5] that the QNL-QC method is not consistent. This obser- 
vation was based on the fact that the authors computed the consistency error in the 
£°°-norm and considered the limit as e — > 0. As a matter of fact, the QNL-QC (and 
even the original energy-based QC method [20]) is consistent in appropriate negative 
norms (see also [TJI1 Lemma 5.11], and [23] for more general discussion). However, 
before we prove this result, we briefly review the notions of truncation error and 
residual, and explain more clearly what we mean by consistency. 

Let y G argmin $ tot (^), then the truncation error (associated with this solution) 
is the linear functional T G ffl*, which measures the extent to which y fails to 
satisfy 

T[u] = D*%(y)[u] = D<S> qc (y)[u] - (f,u). 

Conversely, if y qc G argmin$*° t (^) then the linear functional R G that measures 
the extent to which y qc fails to satisfy ^ is called the residual (of the approximate 
solution y qc ): 

R[u] = D<S> tot (y qc )[u] = D$(y qc )[u} - (/,«). 
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Since y satisfies ^ and y qc satisfies ^ we deduce that 

T = D% c (y) -D$(y), and R = D$(y qc ) - D% c (y qc ). 

Thus, estimating the truncation error will automatically give us a residual estimate 
and vice-versa. Hence, we will call both R and T simply the consistency error. (This 
observation ceases to be valid if coarse-graining is taken into account.) 

Traditionally, we call a numerical method consistent if the truncation error tends 
to zero, in a suitable sense, as the mesh size, or another discretization parameter, 
tends to zero. In the present case, since e is fixed (it is best to think of e to be in of 
the order O(10 -3 ) or O(10~ 4 )), we cannot make such a definition. We will therefore 
use the term rather loosely and simply discuss the order of consistency of the method. 
For example first order consistency would mean that the consistency error can be 
bounded by e and a factor that depends on the smoothness of y. However, even this 
nomenclature is not wholly appropriate as the QNL-QC method has different orders 
in different parts of the domain. This is demonstrated in the following theorem, 
which provides a sharp bound on the consistency error of the QNL-QC method in 
the ^~ 1,p -norms, 1 < p < oo. Even though we will later only use ^ _1 ' 2 estimates, 
the more general case is included here since it requires no additional work. 

Theorem 2 (Consistency in ^ 1,p ). Let y 6 <3f such that min^;^ > 0, then 

\\D$(y) - £>$ q c(y)lk-i* < ^ 2 \\y"\\ im + e*C 3 (W"\\t, { cw + Il3/"ll# (c) ), 
where Ci = Cj(2min^ e c' i = 2,3. 

Proof. The first variations of the atomistic and QNL-QC functionals at y are, respec- 
tively, given by 

N N 

D$(y)[u} = eJ2<f>'(ytH + + + ( 12 ) 



and 



5=1 e=i 



N 



u' e + £ ^ <?%g + y' £+l ) [u' £ + u' £+1 ] (13) 
+ e£§{ 2 W«K + 20'(2<4 +1 K e+1 }. 
The consistency error T := DQ(y) — DQ qc (y) therefore satisfies 

t[u] =e y £W(y'^ + y^M + ^ +1 } - W e K- W^iK+i} ( 14 ) 

Sec 



for all u G *% . From (14) we can immediately obtain a first-order consistency error 
estimate, however, one can improve upon this by taking into account symmetries of 
the interaction law. 

If we collect all coefficients that premultiply a term u'f, then we obtain second- 
order errors in the "interior" of the continuum region and first-order errors in the 
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interface. To see this, we recall the definitions of the interface regions from Section 



2.2, and compute 



(15) 



T[u] = e E {W* + l4fi) " W €+1 )K+i + e E { Wc + *4+i) " W € )K 

+ ^ E {^(^-i+^)+^(^ + ^ + i)-20'(2^)K- 

fec\x' 

Taylor expansions of the first and second terms in curly brackets yield 

fiiv't + y'z+i) - ^( 2 ^+i) = - £ 0"( 2 ^+i)y| + bV'^u)!^'! 2 , and 

where i?^ G conv{^ + y'^ +v 2y^ +1 } and i? 2 ,g e conv{^ + y^ +1 ,2y'^}, and hence 
!</>'"(%) I < C 3 (2min^ c ,^) =: C 3 for i = 1,2. 
To expand the third term in curly brackets in (15), we use (16) twice to obtain 

+ y't+i) + fiiv't-i + y[) - 20'(2j4) 
= eVW*)*? + e 2 \{<p"\^)\y'l\ 2 + <m,e-iM-il 2 }- 



(16) 



Inserting these expansions and bounds into (15), and applying several weighted 
Holder inqualities we can estimate 

\T[u]\ < eJ2{^\v'l\ +^3|i£| a }K + il +sJ2 {eCiWH +e a ^ s |»?| 2 }KI 
£ec\x r 



< {sC 2 \\y"\\ m) + e 2 C 3 \\y'"\\ encxIr) +e?C 3 \W'\\ 2 & {c) W 



□ 



Remark 2. The consistency error has no contribution from the atomistic region, 
which is not surprising since the model is exact in A. In the continuum region, the 
consistency error reflects the second-order accuracy of the Cauchy-Born approxima- 
tion for simple lattices [IT] . Finally, in the interface, the consistency error is only 
of first order. By estimating Hj/'H^m above by 1 1 2/ 7/ 1 1 ^ °° (x) ? we can however obtain 
additional powers of e, 

eC 2 \\y"\\ im < e 1+1 /^ 2 (#J) 1 / p ||y // ||^ ( x), (17) 

but at the cost of a dependency on the size of the interfacial region. Moreover 
it should be stressed that such an estimate might hide the significantly reduced 
accuracy of the QNL-QC method in the interface. 

This loss of accuracy at the interface can potentially become problematic if very 
high accuracy is sought. For example, a P2 finite element discretization of the contin- 
uum region will not be able to increase the accuracy of the computation beyond the 
error introduced in this interface. Thus, a higher order coupling mechanism would 
be a highly desirable goal. 
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A finer analysis of this reduced accuracy for a linearized model problem can be 
found in [5]. □ 

Remark 3. Generalizations of Theorem [2] for a variety of QC methods can be 
found in [23]. There, we also include finite element coarse graining in the analysis. 
We use negative norm consistency error estimates, similar to Theorem |2j to control 
what we call the model error. □ 



Remark 4- In [T9l Lemma 5.11] a negative-norm truncation error estimate with 
respect to a so-called Spijker norm can be found, which, in the language of the present 
work, reads 

\D^(y)[u]-D% c (y)[u]\ < Ce(e- 1 ' 2 \\u'\\ e ^. 

This estimate follows immediately from (17). The reason for choosing this norm is 
that it guarantees uniform closeness of gradients (i.e., e -1 ' 2 ||w'||^| > Hm'I^oo), a fact 
that is highly useful in the nonlinear analysis. Unfortunately it leads to suboptimal 
error estimates. In the analysis of Section [5j we will circumvent this problem by 
using inverse estimates, but still retain the optimal truncation error estimates. □ 



4. Stability 

Having established a fairly sharp consistency error estimate, we turn to the ques- 
tion of stability of the QNL-QC method. Stability estimates are usually easiest to 
establish in a Hilbert space, particularly, when the operator under investigation is 
coercive. Coercivity of the QNL-QC hessian evaluated at the reference state y = Fx 
was established in [5j [19], and it was shown in [9] that these estimates sharply reflect 
the stability of the full atomistic model. In the following sections we will generalize 
these results to various situations that admit large deformations. 

In one dimension, it is still relatively straightforward to derive stability bounds 
in the space < % 1 ' p . For example, several stability results for atomistic and quasi- 
continuum models in the space ^f 1,QC , which are particularly useful for a nonlinear 
analysis, were derived in [22]- However, since such results would be very difficult 
(if not impossible) to obtain in 2D/3D, we will attempt in the next sections to use 
only ^ 1,2 -stability results in our analysis and still obtain similarly strong nonlinear 
results. 



4.1. Preliminary remarks. Linear stability for minimization problems is normally 
associated with the coercivity constant (or smallest eigenvalue) in an appropriate 
function space. Thus, to relate the stability of the atomistic model and that of the 
QNL-QC, we would like to prove a result of the type 

c(y) := inf D 2 $(y)[u,u] « inf D 2 & qc (y)[u,u] =: c qc (y). (18) 

11^11^1,2=1 1,2=1 

If the approximate hessian -D 2 $ qc could be obtained by perturbing the coefficients 
of D 2 Q then one could obtain such a result from an error estimate on the hessians, 
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that is, from a bound on \\D 2 §(y) — D 2 ^ qc (y)\\ L ^i,2 ^-i,2). However, if we compute 
the two hessians explicitly, 

N N 

D 2 $(y)[u,u] = e^(j)"{y^)\u^\ 2 + e <p"(y'^ + y' i+1 )\u'^ + u^\ 2 , and 

N 

D 2 % c (y)[u,u]=eJ2<P"(ydK\ 2 + e£W* + *4+i)K + 

+ £ £i{0"(224)|24| 2 + 0"(2^ +1 )|2<4 +1 | 2 }, 

we observe that in the continuum region D 2 Q(y)[u, u] contains non-zero coefficients 
for the mixed terms u^u^ +1 which are replaced by diagonal terms in the QNL-QC 
Hessian. This shows that it is impossible to obtain a useful estimate on the difference 
of the Hessians. Instead, we must aim to prove ( 18 ) directly. 

The crucial observation which is the basis of the stability analysis in this section 
is that the non-local hessian terms \u'^ + ?4 + i| 2 can be rewritten in terms of the local 
quantities \u'^\ 2 and |w^ +1 | 2 and a strain-gradient correction, 

\u^ + u^ +1 \ = 2\uA +2\u£ +1 \ —e \u^\ . (19) 
Using this simple formula it is straightforward to rewrite the Hessians in the form 

N N 

D 2 <S>{y)[u,u] = s^A ( \u' ( \ 2 + e^e 2 B ( \ul\ 2 , and 

D 2 % c {y) [it, u\ = e E^KI 2 + e E ^eKI 2 . 
£=1 

where 

A<: = (f>"(y^ + 20"(x4_i + <4) + 20" + 

f 20"(^-i + vd + 20" + £ e .A', 

20"(^ + ^ +1 ) + 20"(2^), 

I VW € ), eeC\l', 

^ = -r(i4+i4fi)- 

The proof of (20) requires purely algebraic manipulations. 

We observe that, apart from estimating the effect of replacing A^ by A%, which is 
the purpose of the next lemma, our main challenge will be to understand the effect 
of dropping the strain gradient term in the continuum region. 

Lemma 3. Let y G & such that min^ y^ > 0, then 

— ^Wt 00 — e/ LC 3 \\y"\\e°°(c), 

where C* 3 = 6*3(2 min^ e c y'g). 
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Proof. For £ G A', A% — A^. For (gC'\ X', we have 

A e - ^ = 2(0"(s4_ 1 + $ - 0"(2^)) + 2(0"(^ + - 0»(2$) 

where G conv{yl 2 + 2?/^} and $2,£ £ conv{y^ + y'^ +1 , 2y'^}. Hence, we obtain 

| ^ - As\ < 2eC 3 \y'l_ 1 \ + 2eC 3 \$\ < 4eC 3 \\y"\\ c . 

Performing a similar calculation for the interfaces, and sacrificing a factor two there, 
we obtain the stated result. □ 



Remark 5. If G C 4 ((0,+oo]) then the result of Lemma [3] can be improved. 
Using similar calculations as in the proof of Proposition [2j one obtains 

HA: - MU™ < max {e2C 3 \\y"\\ i oc i j ) ,e 2 C A (\\y" , \\^ {c > ) + \\y"\\%°(C)))- 
Even though this result is clearly sharper than Lemma |3j it carries only limited 
practical value due to the fact that the maximum in ||y"||£°°(e) wm typically be 
attained at the interface. □ 

Before formulating the first stability result, we briefly discuss the nature of the 
coefficients B^. In view of the discussion in Section 2.3 we have assumed throughout 
that y'^ > r*/2 for all £ G Z and consequently, 4>"(y^ + v'm) < for all £, that is, we 
obtain the following result. 

Lemma 4. If y'^ > r*/2 for all (6Z, then 

B^>0 V£ G Z. (21) 



4.2. An illustrative example. Before we turn to the stability estimates, we briefly 
discuss an example that further stresses the difference between atomistic and QC 
hessians. We consider the case y = Fx and .4. = 0, that is, the QNL-QC method 
reduces to the local QC (or simply, Cauchy-Born) method. In that case, we have 

D 2 $(Fx)[u,u] = A\\u'\\ 2 t 2 +e 2 B\\u\\ e z and D 2 $ qc (Fx){u, u] = A\\u'\\%, 

where A = <j>"{F) + 40"(2F) and B = -<j>"(2F). 

In this case, the ^ ^-spectrum of D 2 § qc (Fx) contains only the eigenvalue A and 
every vector u G % is an eigenvector. On the other hand, the ^ ^-spectrum of 
D 2 ^{Fx) is given by 

%i,2(D 2 $(Fx)) = {A + BfjLj : j = 1, . . . ,N - l}, where /ij = 4 sin 2 (fjW); 

see [9l Sec. 3.1] and P, Lemma 1] for similar analyses. 

Thus, we see that the eigenvalues of the low eigenmodes are well approximated by 
the local QC stability constant A, and in particular, 

c(Fx) =A + 4B sin 2 (vre) = A + 0(e 2 ) = c qc (Fx) + 0(e 2 ). 

On the other hand, the eigenvalues corresponding to the high frequency modes are 
significantly larger, in other words, high frequency perturbations are considerably 
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more expensive in the atomistic model than in the continuum model. This discrep- 
ancy is not too surprising since the basic assumption of the Cauchy-Born model is 
the absence of high-frequence modes in the deformation. 

4.3. An a priori stability result. A common approach to establish the stability of 
a QC method is to assume that next-nearest neighbour interactions are dominated by 
nearest-neighbour interactions jH [151 EE 122] . This leads to particularly simple results 
that can even be made fairly sharp as the numerical experiments in [22] demonstrate. 
If one wishes to include more than the simplest examples in such stability results, 
then one should specify classes of atomistic solutions that are of interest and prove, for 
each such class, that the QC hessian is stable when evaluated at those deformations. 
In [22], two main types of stable solutions were identified for a variant of the model 
problem (J3|: (i) elastic deformations and (ii) deformations with a single crack. Here 
we will focus only on elastic solutions only, for which we obtain the following result. 

Proposition 5 (Stability of elastic states). Let y e W , such that y'^ > r*/2 for 
all (gZ, then 

c(y) > min A%, and c qc (y) > min A%, (22) 



where c(y),c qc (y) are defined in (18) and A^A^ are defined in (20) 



Proof. Both bounds in (22) follow immediately from (20) and from the fact that 



> for all ( 6 Z. □ 

Remark 6. Proposition [5] provides a particularly simple example of atomistic 
configurations for which the QNL-QC method is stable. Due to the representation 



of next-nearest neighbour interactions developed in Section |4.1[ these results are 
in fact significantly sharper than those found in [H [THl 122], and will be sufficient 
whenever elastic effects dominate. Since there is no regularity assumption on y, 
simple singularities are included as well. □ 

Remark 7. It must be said at this point that the stated goal, namely a proof 
of c(y) ~ c qc (y), or at least c qc (y) > c(y), up to a controllable error, has not been 
obtained here (however, a reverse inequality will be established in the next section). 
It can be easily seen that such a result would be related to the regularity of eigen- 
functions in the continuum region. Namely, let u be an eigenfunction for D 2 § qc (y), 
associated with the eigenvalue c qc (y), then 

c(y) = min D 2 &(y){u, u] 

ll M 'll/2=l 

< D 2 ®{y)[u,u] 

< c qc (y) + \D 2 $(y)[u,u] - D 2 <S> qc (y)[u,u]\. 



Using (20) and Lemma |3j we obtain the estimate 

\D 2 ^y)[u,u]-D 2 % c (y)[n,u]\ < iC 3 e\\y"\\^ {Q +e 3 J2 B H 



u"\ 2 . 
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Thus, if one could establish an 0(1) bound on ||w"||^(c)5 f° r example in terms of the 
smoothness of y in the continuum region, then one would obtain 

c qc (y)>c(y)-0(e). 

However, such a bound is far from trivial to obtain, even in one dimension. □ 

4.4. A posteriori stability. If we assume that y = y qc is a local minimizer of the 
QNL-QC approximation, then we can check the condition whether D 2 Q qc (y qc ) > 
a posteriori by solving a (generalized) eigenvalue problem. The question we then 
ought to ask is whether this stability carries over to the full atomistic model, that is, 
whether also D 2 §(y qc ) > 0. If this were not the case then it cannot be guaranteed 
that a local minimizer in the QNL-QC model corresponds to a local minimizer in the 
atomistic model. However, for our model problem the answer is surprisingly simple 
and is given in the following result. 

Theorem 6 (A Posteriori Stability). Let y G <3f such that mm^zy^ > 0, and 
such that 

Wj + 4i)<0 V£gC, 

then 

c(y) > c qc (y) - eAC 3 \\y" Weo,^, 

where C 3 = C 3 (2 min feC , 

Proof. From Lemma [3] and the assumption that (f>"(y'^ + y'^+i) < 0, or equivalently, 
Bj: > for all £ G C we have, for all u G 

$"(y)[u,u] > ^MWM - e4C 3 \\y"\U~ {e) \\u'\\ Ps . 

The stated result follows by taking the infimum over tiff' with ||u||^i,2 = 1. □ 



Remark 8. The analysis in [13] shows that Theorem |6j is false in higher dimenions. 
Namely, it is possible that a QC model (for example, a pure Cauchy-Born model) 
is stable at a uniform deformation Fx while the atomistic model is unstable. Thus, 
at the very least, one needs to add an additional assumption that the deformation 
gradients in the continuum region belong to an atomistic stability regime. □ 

4.5. A remark on defects. Suppose now that a deformation y 6 f has a single 
'crack', that is, there exists a single index £ such that y'~ = O(N). Assume, for 
example, that F > and that y G W is given by 

" — / 1, if ^ 7^ ^ / 2 o\ 

y ^ ~ \ F + (N - l)(F - 1), otherwise. 1 ' 



If iV is sufficiently large, then < and hence Proposition (22) gives no useful 



information. This is no accident as the following result demonstrates. 
Proposition 7. Suppose that there exists a cut-off radius r c > r* such that = 



in [r c , +oo). Let y G & be defined by (23), then 



<Ae and c qc (y) < Ae, where A = 0"(1) + 40"(2). (24) 
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Proof. There exists a unique test function such that 

(iV-l) 1 ^, £ = 1 



**~ 1 -(iv-ir 1 / 2 , 

and we note that = 1. 

Next, we observe that u'^ = for all £ except for £ = £, £ — 1. However, since = 
in [r c , +oo), it follows that B^_ x = = for iV sufficiently large, and hence we can 
simply ignore the strain gradient term. Thus, testing D 2 Q(y) with u, we obtain 

N 

D 2 <t>(y)[u,u] = eJ^MN ~ I)' 1 - 

Setting A = 0"(1) + 40"(2), we have ^ = A for all £ except £ = | - + 1. 

For £ = £, we have A| < provided that iV is sufficiently large, while, for £ = £ ± 1, 
it follows from the concavity of in (r*, +oo) that < A. Hence, we obtain 

JV 

-D 2 $(y)[-u, tt] = e^A ? (iV- I)" 1 < ei. 

Since c(y) > c qc (y) for this example, the same result also holds for the QC stability 
constant. □ 

Since deformations with cracks or other defects can surely be stable equilibria, 
this result only shows that they cannot be analyzed in the ^f 1,2 space. Different 
topologies must be chosen to analyze the stability of defects. The analysis in [22] 
suggests that the correct norm for a ID example with a single crack might be 

1/2 



' |2 



however, we will not pursue this direction further in the present work. We note, 
however, that changing the norm inside the atomistic region does not change the 
residual estimate. To this this, one verifies that, in the last line of the proof of 



Theorem 



the term ||w'|| £ p' may be replaced by ||w'||^/, c ,j. 

5. Existence and Convergence 



In this section, we will combine the consistency and stability analysis of Sections 
[3] and [4] to prove a priori and a posteriori existence results and error estimates. 

5.1. An a priori existence result for elastic deformations. In this section, we 
extend the main results in [5, 19] to the context of large deformations and solutions 
that are possibly non-smooth in the atomistic region. 

We begin by assuming the existence of a strongly stable equilibrium y of the full 
atomistic model. If we assume that this deformation is "sufficiently smooth" in the 
continuum region then Theorem [2] shows that the truncation error is small, and hence 
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y can be considered an approximate solution of the QNL-QC equilibrium equations 
If ?/ is also stable in the QNL-QC model, then we can employ the inverse function 
theorem, Lemma [TJ to prove the existence of an exact solution y qc of the QNL-QC 
equilibrium equations in a small neighbourhood of y. This idea is made rigorous in 
the following result, which is best read as follows: 

"If y is a strongly stable local minimizer o/$ tot ; which is smooth in the continuum 
region, then there exists a solution y qc of the QNL-QC method, which is a good 
approximation to y." 

Theorem 8. Let y G argmin$ tot ; and assume that > r*/2, and that 

min A £ =: A > 0. (25) 

Then there exist constants 61,62 that depend only on min^y^ and on A, such that, if 

£ \\y"\U°°(C) < 61, and (26) 

e 1/2 \\y"\Um +e 3/2 (\\y'"\k(c\T ) + \\y% {c) ) < (27) 

then there exists a (locally unique) local minimizer y qc of in *3/ such that 

\\(y-y qc )'\k < i{eC 2 \\y"\\ m) +e 2 C 3 {\\y>"\\ el{c ,\ x ,) + \\y"\\l ei c))}- 
where Ci = Cj(2min^ eC / y'^), i = 2,3. 

Proof. To put the statement into the context of Lemma [T] we define X = % 1,2 , 
Y = ^T" 1 ' 2 , &{w) = D&°*(y + w), and x = 0. The set A is given by 

A = {w G : ^ + w'z > 0, £ G Z}. 

Since (ft G C 3 ((0, +oo]), it is clear that & is continuously (Frechet) differentiable, 
with D3?(0) = D 2 <f> qc (y). 

1. Consistency: From Theorem [2] we obtain that 

||^(0) \\y = \\DQ™{y)y-i,2 = \\D% c (y) - D$(j/)||*-i, a < 77, where 
T) = eC 2 \\y"\\ em +e 2 C 3 (\\y'"\U2 ( c>\i>) + h"\\%(c))- 

2. Stability: The assumption that y'^ > r^/2 implies that B £ > for all £. 
Using this fact, the "elasticity assumption" (25), as well as Lemma [3j we obtain 

N 



D 2 <S> qc {y)[u,u] > eJ2M u \\ 



€=1 

N 

> eJ2 A aK\ 2 - £ 4C 3 \\y"\\^ iC )\\u'\\ 2 £m 

> (A-eAC 3 \\y"\U^ C ))\\u%- 



Hence, we obtain 



c q c(y) >A-e4C 3 \\j/'\ 
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and consequently, if A — e4C 3 \\u"\\^(c) > 0, then 

P^Or 1 !^-!,^) < (A-eAC 3 \\y"\\^ (c) y 1 . 
Setting 5 l = A/{8C 3 ), @ becomes 

e4C 3 \\y"\\ tOB( c) < \A, (2? 

and hence we obtain that 



D^(0)- 1 \\ m -i,2 ! ^,2 ) <a:=2/A. 



3. Lipshitz bound: Next, we compute a Lipschitz bound for = D 2 Q 



qc • 



To this end, we first need to establish a lower bound on y'^ + w'^ for all w satisfying 
< 2r]a. Since the ^ 1,2 -topology cannot provide this bound directly, we need 
to resort to an inverse inequality: for all w G such that ||iy||^i,2 < 2t]<j, we have 

IKHi- < e-^Wv/Wz < M[(eV*\\y"\\ m) + e^{\\y'"\\ Pe{c , xv) + \\y% (e) )), (29) 

where M[ is a constant that depends only on A and on C 2 ,C 3 . Hence, if 

e 1/2 \\y"\Ui(i) + e^WkicW) + \\y%(c)) < w { ™*l/ v (30) 

then the desired bound 

y't + w'z > Sy't 

holds. 

It is now a simple exercise to show that, for any two displacements wi,W2 G & , 
such that ||wj||^i,2 < 2rja, we have the Lipschitz bound 

\D 2 <$> qc (y+w 1 )[u,v]-D 2 <$> <ic (y+W2)[u,v}\ < L'\\(w 1 -w 2 )'\\^\\u'\\ e 2\\v'\\ l 2 W,v G <2r, 

where V is a constant that depends only C 3 (Smm^y'^). Applying the same inverse 
inequality as in (29) we obtain that, 

\\DJ?(wi) - DJ?(w2)\\L(wh2,<&-i,2) < e~ 1/2 L'\\wi - w 2 \\^i,2. 

We set L = s~ l l 2 L' . (At this point we could choose 5 in order to optimize the size of 
L against the requirement (30). Since we are only interested in a qualitative result, 
we will ignore this possibility.) 

4- Conclusion: To conclude the proof, we only need to check the condition 
under which 2La 2 r] < 1. Inserting the bounds for the various terms, assuming that 
(30) and (28) hold, we obtain that 

2Lo\ < M' 2 {e^\\y"\\ Em + E^(\\y"%> { o\v) + ||/||| (c) )} 

where M' 2 depends only on A, on C 2 , C 3 , and on L'. 
Hence, we conclude that, if (26) and (27) hold with 

<Ji = A/(8C 3 ) and 5 2 = min min^, l/M' 2 }, 

then there exists w qc G % such that F(w qc ) = 0, or equivalently D^„°^(y + w qc ) = 0, 
and 

||w qc |ki,2 < 2r]cr < jT]. 
Taking into account the bound for r\ from step 1, we obtain the stated error estimate. 
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Local minimality of y qc = y + w qc follows from ( 26 ) which implies that 

Cqc(?/) > | A 

and the local Lipschitz bound, which guarantees 

c qc (?/qc) > c qc (y) - L(2rja) > c qc (y) - \ = 0. 
Thus D 2 Q qc (y qc ) is positive definite and therefore y qc is a local minimizer. □ 



Remark 9. As indicated above, Theorem 26 is a variation or generalization of [5j 
Thm. 4.1] and [19], Thm. 5.13], and hence we should briefly discuss the differences 
between these results. 

Theorem 4.1 in [5] treats a linearized case, but is otherwise fairly sharp. For 
example, the stability condition derived in [5J Lemma 2.3] was later shown to be exact 
[§]. The convergence rate in terms of the atomic spaces e (or h in p)j) is optimal, 
however, the dependence of the error on the smoothness of the exact solution is 



much stronger than in Theorem [26J This is a consequence of deriving the trunction 
error estimate in an £ p -type norm, and comparing the QC solution with a continuum 
solution instead of the exact atomistic solution. 

Although a very different terminology is used, some aspects of the proof of [191 
Thm. 5.13] are quite similar to the analysis presented here. However, in [19] the 
consistency analysis is based on a higher order expansion of the atomistic solution 
in terms of the Cauchy-Born model. As a result, the analysis proves the existence 
of atomistic and QC solutions only for "sufficiently small and sufficiently smooth" 
applied forces, in a small neighbourhood of the reference state Fx. Moreover, the 
use of a different norm for measuring the truncation error [TH] Eq. (5.27)] leads to 
error estimates that are optimal in e only in the ^ ,1 '°°-norm, and it is unclear from 
the presentation how the error depends on the smoothness of the atomistic solution. 

The work in [19] is more general than the present paper in that it analyzes the 
more general geometrically consistent QC scheme proposed in [10] . □ 



5.2. An a posteriori existence result. Theorem [8] is an a priori result, that is, it 
guarantees the existence and accuracy of a QNL-QC solution under certain assump- 
tion on a given atomistic solution. Results of this type can guarantee that certain 
classes of atomistic solutions can be reliably approximated by the QNL-QC method. 
In practise, however, an a posteriori result of the same type may be even more in- 
teresting. Given a computed QNL-QC solution y qc , we may ask whether an (exact) 
atomistic solution y exists that y qc is an approximation to. Questions of this kind are 
investigated in the literature on numerical enclosure methods (see [21] for a recent 
review, or [21] where this concept is called a posteriori existence) and leads to the 
following result, which is best read as follows: 
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"If y qc is a strongly stable QNL-QC solution, which is smooth in the continuum 
region, then there exists an exact solution y of the full atomistic model, for which y qc 
is a good approximation. " 

Theorem 9. Let y qc £ arguiin^*^, such that c qc (y qc ) > and 

0"((2/qc)f +(2/qc)e + l) <0 V^£C. 

Then there exist constants 5i,5 2 that depend only on mm^(y qc )^ and on c qc (y qc ) such 
that, if 

£ \\y'L\U°°(C) < 5i, and (31) 



e 



y' qc \Vm + £ 3/2 (llj/qclUi(c'\xo + lbqcll|(o) < ^ (32) 

then there exists a strongly stable local minimizer y of $ tot in IV such that 

||, _ yll < £( ^2||< C ||^(J) +g 2 C , 3 (||<cll^(C'\JQ + lbqcll|(C)) 

Cqc{y qc ) - £4:G 3 \\y h \\ e ^ {C ) 
where Q = d(2 min ?eC /(?/ qc )^), i = 2,3. 

Proof. The proof follows along very similar lines as the proof of Theorem [8j and we 
will therefore only focus on those parts of the argument that change. 

To put the statement into the context of Lemmaflj we define X = Y = ^~ 1,2 , 
JP(w) = D& tot (y qc + w), and Xo = 0. The set A is again given by 

A = {w £ <& : (y^ + > 0, £ £ Z}. 

Since £ C 3 (0, +oo), it follows that & is continuously Frechet differentiable, with 
D&(0) = D 2 <S>(y qc ). 

1. Consistency: From Theorem [2j we obtain 

||^(0) \\ Y = ||£>$ t0 % qc ) ||*-i. a = ||D$(y qc ) - D% c (y qc )y^, 2 < rj, where 

V = £ C 2 \\y qc \\p e( x) + £ 2 C 3 {\\y qc \Ui(c>\i>) + \\y qc \\%(c))- 

2. Stability: Using Theorem [6j we obtain 

c{y qc ) > c qc {y qc ) - e4C 3 \\y''\\e°o( C) =: 1/a, 



where C 3 = C 3 {2 min^c'G/qc)^)- Hence, if we require that (31) holds with 8\ = 
Cqc(y q c)/(8C , 3), then we obtain 

\\D^(0)' l \\ L ^-i,2 t ^i,2) < l/c(y qc ) <a< 2/c qc (y qc ). 

3. Lipschitz bound: The proof of a Lipschitz bound for is a verbatim 
repetition of step 3 in the proof of Theorem [8j For some fixed < 5 < 1, we 
obtain that, if ||^||^| < 2rja, i = 1,2, then 

\\D&(wi) - D^(w 2 )\\l{w^i,w-^) < e^^L'Wwx - w 2 \\^i,2, 
where L' depends only on C 3 (5min^(j/ qc )^). The actual Lipschitz constant is again 
set to L = e~ 1/2 L'. 
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4- Conclusion: The conclusion of the proof follows precisely step 4 in the proof 
of Theorem [8j The different constant in the error estimate is due to the sharper 
definition of a. □ 

Remark 10. Theorem [9] is the first rigorous a posteriori error estimate of this 
type for the QNL-QC method. However, the basic idea behind the result, that is, to 
use a posteriori information about the approximate solution to deduce the existence 
of an exact solution, is not new. A similar approach to the one used in the proof of 
Theorem [9] can be found in the work of Plum [21] on numerical enclosure methods 
for semilinear differential equations. A detailed analysis of these ideas, in the context 
of a posteriori error control for finite element methods was given in [21]. The first 
application of this idea in an analysis of a QC method was presented in [221 Thm. 
5.1]. The latter result differs from Theorem [9] not only in the type of QC method 
analyzed, but also in the choice < W 1, °° for the function space setting. While that choice 
gives uniform neighbourhoods and thus asymptotically sharper existence conditions, 
a sharp stability result such as Theorem [6] would be difficult to prove. Moreover, a 
generalization of [221 Thm. 5.1] to higher dimensions seems virtually impossible. □ 

Conclusion 

In order to outline a new framework for the analysis of QC methods, we have chosen 
a particularly simple setting, admitting only second-neighbour pair interactions, and 
disregarding finite element coarse graining of the continuum region. Extensions of 
this work to more complex interaction models and to coarse-grained models are in 
progress [231 HI"] . The most interesting extension, however, will be to the application- 
relevant two- and three-dimensional setting. While great care was taken in this paper 
that a clear path for such a generalization is available, many interesting questions 
will need to be carefully considered in order to establish the consistency and (linear) 
stability of the QC method. Even the relatively basic question of consistency is still 
open for a general finite range interaction model, however, the work in [TU] provides 
a promising starting point. No results on the stability of the QC method in 2D/3D 
can be found in the literature at present. 
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